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SUMMARY 

A three-dimensional, geometrically nonlinear two-node Timoshenko beam element based 
on the Total Lagrangian description is derived. The element behavior is assumed to be 
linear elastic, but no restrictions are placed on magnitude of finite rotations. The 
resulting element has twelve degrees of freedom: six translational components and 
six rotational-vector components. The formulation uses the Green-Lagrange strains 
and second Piola-Kirchhoff stresses as energy-conjugate variables and accounts for for 
bending-stretching and bending-torsional coupling effects without special provisions. 
The core-congruential formulation (CCF) is used to derived the discrete equations in a 
staged manner. Core equations involving the internal force vector and tangent stiffness 
matrix are developed at the particle level. A sequence of matrix transformations carries 
these equations to beam cross-sections and finally to the element nodal degrees of free- 
dom. The choice of finite rotation measure is made in the next-to-last transformation 
stage, and the choice of over-the-element interpolation in the last one. The tangent 
stiffness matrix is found to retain symmetry if the rotational vector is chosen to mea- 
sure finite rotations. An extensive set of numerical examples are presented to test and 
validate the present element. 
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1. INTRODUCTION 

The computer-based geometrically-nonlinear analysis of flexible three-dimensional structures 
has attracted considerable interest in recent years. In the aerospace field, part of this atten- 
tion comes from establishing challenging and ambitious goals for space research, such as the 
space station, space-based antennas for improved communications, space-based telescopes, so- 
lar arrays, and a large variety of similar devices which have been given detailed consideration 
in various proposals for future space developments. In mechanical^ engineering, interest has 
emerged from the competition to obtain reliable, accurate and inexpensive manufacturing pro- 
cedures, especially in the development of the new generation of robots capable of performing 
high precision tasks such as microspot welding and assembly of orbiting structures. Further 
interest comes from analysis and design of high-performance aircraft, helicopter blades and 
turbomachinery. These applications have motivated the development of more physically re- 
alistic computational models of large flexible structures that exhibit pronounced geometric 
nonlinearities. 

Three kinematic descriptions have been used in geometrically nonlinear finite element 
analysis: Total Lagrangian, Updated Lagrangian and corotational. The present work follows 
the Total Lagrangian (TL) description, but in an unconventional variant that constructs 
the nonlinear finite element equations in a staged fashion. This variant is called the core- 
congruential formulation and identified by the acronym CCF in the sequel. An account 
of this methodology is presented in a recent review paper by Felippa and Crivelli. 1 This 
review concludes that a key advantage of the CCF for constructing TL elements is that it 
helps establishing consistency by avoiding the premature introduction of drastic kinematic 
approximations. 

The main ideas behind the CCF can be traced to a 1973 paper by Rajasekaran and 
Murray 2 , who examined critically the pioneer work on the Total Lagrangian description 
by Mallett and Marcal 3 . The 1974 discussion of Rajasekaran-Murray’s paper by Felippa 4 
established general expressions for the finite element equations that appeared at various 
variational levels. Further historical details are given in the review by Felippa and Crivelli. 1 

This paper uses the CCF to derive the finite element equations of a TL three-dimensional 
Timoshenko beam element that can undergo arbitrarily large rotations. First, we derive the 
governing differential equations encountered in the geometrically nonlinear static structural 
analysis of three-dimensional beams. Next, the finite element counterparts are obtained by 
discretizing the physical degrees of freedom. Our main assumption is that the beam behaves 
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as a linear hyperelastic isotropic medium, which allows us to write its internal energy as 
a quadratic function of the finite strains. We obtain the equilibrium equations from the 
stationarity condition on the first variation of the total energy. Similarly, we obtain the rate 
or incremental equations from the second variation of the total energy. 

The CCF derivation of the governing equations of motion proceeds through two phases: 
a core phase followed by a transformation phase. In the initial phase core energy, residual 
and tangent stiffness matrices as well as internal force vectors, are obtained independently of 
any specific choice used to represent or parametrize the motion. These matrices and vectors 
pertain to individual particles. They do not depend on discretization decisions, such as 
element geometry, shape functions and selection of nodal degrees of freedom. To emphasize 
this independence, the term core was coined. In the transformation phase, these core forms 
are gradually specialized to particular element instances. This specialization is achieved by 
the application of one or more transformation stages that progressively “bind” particles into 
lines, areas or volumes through kinematic constraints, and eventually link the element domain 
to the nodal degrees of freedom. The choice of specific parametrizations for finite rotations 
may be deferred to latter stages. 

What are the differences between the CCF and the more conventional Total Lagrangian 
formulation of nonlinear finite elements? If kinematic exactness is maintained throughout, the 
final discrete equations are identical. But in geometrically nonlinear analysis approximations 
of various kinds are common, especially in structural elements with rotational degrees of 
freedom such as beams, plates and shells. In the conventional, one-shot formulation it is 
often difficult to assess a priori the effect of seemingly innocuous approximations “thrown 
into the pot,” and a posteriori exhaustive testing of complex situations becomes virtually 
impossible. Sample: how does the neglect of higher order terms in the axial deformation of 
a spinning beam affects torsional buckling? 

The staged approach taken in the CCF permits a better control over such assumptions. 
The core equations Eire physically transparent, clearly displaying the effect of material be- 
havior, displacement gradients and prestresses. In the ensuing transformation sequence the 
origin of each term can be exactly traced, and on that basis informed decisions on retention 
or dropping made. 

Another important advantage of the staged approach is the precise identification (and 
avoidance) of kinematic choices that lead to unsymmetries in the tangent stiffness. In beams, 
plates and shell elements such a symmetry loss is linked to the choice of the finite-rotation 
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parametrization (Euler angles, Euler parameters, rotational vector, direction cosines, etc). 
This decision can in fact be postponed to the final stages, and changed if necessary with- 
out affecting “kernel” forms derived in previous stages. This nesting has obvious beneficial 
influence on element programming modularity. 

2. PREVIOUS WORK IN NONLINEAR SPACE STRUCTURES 

Early work in the area of orbiting and free-flying structures dealt with the problem of rigid 
spacecrafts with flexible appendages attached to their core. 5 In this case, the motion of the 
system was obtained by superposing a given number of linear elastic modes to the overall 
motion obtained considering the structure as made of interconnected rigid bodies. This 
procedure has been called the hybrid coordinate method. 6 The flexible motion — assumed 
small — is then described with respect to frames attached to the underlying rigid core motion. 

When this procedure is to be applied to structures with distributed flexibility — those not 
having a distinct rigid core — a question that immediately arises is how to choose the reference 
frame. One idea is to define a floating or unattached frame that is optimum in some sense. 
Two frames were proposed by De Veubeke, 7 one that minimizes the relative kinetic energy 
and another that minimizes the deformation energy. The first choice is shown to correspond 
to Tisserand’s conditions of zero relative momentum and angular momentum, i.e., the rigid 
body modes are found to be fixed in this frame. However, this choice introduces some practical 
difficulties, especially in the case where there are lumped sources of kinetic energy, such as 
rotating masses or gyros. Further work has also been done to include the effects of spinning 
rotors. 8 The hybrid coordinate method has the advantage that the equations of motion are 
represented in a form similar to rigid body equations. 9 This type of moving frame has also been 
called mean axis system 10 and used to implement finite elements representations. Extensions 
of the mean axis formulation to flexible multibody system dynamics can also be found. 11 De 
Veubeke 7 shows that the use of the minimum deformation energy criterion allows the relative 
displacement to be exactly represented by an expansion in natural elastic vibration modes 
and leads to a simpler implementation. In any case, introducing a floating frame requires 
constraint equations to be added, because they require the definition of additional variables 
that cannot be obtained directly from the dynamics of the system. 

A more subtle problem associated with the floating frame formulation is found in con- 
nection with spinning free-free beams. 12 In this case, nonlinear effects produce a geometric 
stiffening due to the spin-induced longitudinal stretch. The resulting axial force, which cannot 
be considered infinitesimal, affects the beam bending stiffness, showing that the uncoupling 
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between rigid body modes and elastic modes no longer holds. Several attempts have been 
made to include these and other effects by higher order corrections to the theory . 13 This 
leads to cumbersome expressions for the matrices entering the governing equations of motion, 
without apparent advantages with respect to a full nonlinear theory. 

With increasing interest in control-structure interaction, the floating frame approach has 
gained new attention, especially when the flexible component is attached to a large rigid mass 
and there is a hierarchical control system that keeps the elastic deformation small. A typical 
example is that of a flexible beam or antenna attached to the space shuttle . 14 In this case, 
the shuttle can be regarded as a rigid body to which the reference frame is attached, while 
the flexible part is discretized using finite elements. The relative equations of motion of the 
flexible part axe linear whereas the nonlinearity comes from the coupling with the rigid body 
motion. Assuming that the inertia of the flexible part is small compared to the inertia of 
the shuttle, the flexible motion can be regarded as a perturbation to the rigid body motion. 
This perturbation technique allows the analyst to define a rigid-body maneuvering strategy 
independently of the elastic behavior. The linearity of the elastic component is required to 
construct an optimal feedback control scheme for vibration arrest. This methodology has 
been used to model the SCOLE experiment . 15 

When performing large-rotation dynamics analysis using the floating frame approach, it is 
important to note that the coupling with the rigid body motion must come through the inertia 
components, because the deformation components have been intentionally uncoupled from 
it. Thus, the inherent nonlinearity of the problem carries over to the inertia terms, leading 
to fully populated nonlinear mass matrices that ruin the sparsity property of conventional 
finite element analysis. For large systems, this is computationally inefficient and prohibitive in 
terms of storage, forcing analysts who use this approach to look for reduction methods. Thus, 
most of the programs based on this approach use linear modes of the free-free structure in 
the undeformed configuration to condense the problem to a few degrees of freedom. Another 
approach is to use a fixed frame for the inertial terms only . 16,17 

Although the single frame approach has been used extensively in spacecraft, its appli- 
cations are limited to small elastic deformations and thus mainly confined to the modeling 
of free-vibration dynamics. When the primary concern is large deflections, as in the case of 
stability and/or postbuckling analyses, the relative rotations between structural components 
are no longer infinitesimal. It is also desirable to preserve the structure of the finite element 
equations for problems of dynamic instability. This has motivated the development of the 
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corotational description. This description retains multiple reference frames relative to which 
the elastic deformation of portions of the structure are described. 18 In its combination with 
finite elements each one has a co-moving attached Cartesian frame. The motion of this rigid 
frame is used to decompose the element displacements into rigid body and deformational 
components. Because the latter are assumed small in each element, small relative deforma- 
tion measures may be used. This model is intimately related to the discretization process, 
i.e., the finite element discretization is done before the establishment of the equations of 
motion or the definition of the variational principle. This procedure has been extended to 
dynamics analysis of space frames 19 but only limited to explicit integration schemes, which 
do not require the explicit computation of the stiffness matrix. 

Because of the small relative deformation assumptions, the expressions of the finite element 
matrices in the corotational frame can be those corresponding to a linear finite element model, 
optionally corrected by geometric stiffness effects. An interesting question that may be raised 
is: Is it possible to obtain a set of external transformations that project these matrices into 
the global frame? This will have the advantage that existing finite elements can be taken as a 
“core” component which can be transformed to the global equations by appropriate external 
manipulations. The answer given by Rankin and coworkers is partially positive. 20,21 This was 
done by enforcing rotational invariance of the internal force, which in turns translates into the 
satisfaction of rotational equilibrium. This technique relies on the use of a projector operator 
which removes the rigid rotations. However, the kinematics properties of the corotated frame 
still depend on a subset of element properties such as dimensionality and the number of nodal 
points. 

Another way to achieve the projection goal is to use a finite strain theory from the out- 
set. In this case, the effect of large rotations is automatically taken into account. Simo 
and coworkers 22 ’ 23 and Cardona 24 have exploited the first Piola-Kirchhoff (PK1) stress, for 
which the conjugate strain is simply the deformation gradient. This leads to a relatively 
straightforward formulation of the discrete equilibrium equations, from which an incremental 
solution procedure is obtained. Downer Park and Chiou 25 have constructed a corotational 
formulation based on Cauchy (true) stress increments and applied to the dynamic analysis of 
spinning beams, with emphasis on energy and momentum conserving algorithms. 

The present work differs from previous ones in the following respects: 

1. The Total Lagrangian (TL) description is used for 3D Timoshenko beam elements in 

conjunction with the second Piola-Kirchhoff (PK2) stress and the Green-Lagrange (GL) 
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strain. A symmetric tangent stiffness is obtained for a particular choice of the finite- 
rotation measure. 

2. No kinematic restrictions axe placed on the overall rotations. Only a mild restriction 
applies at the element level: the relative rotations within an individual beam element 
should not exceed 360°. 

3. The CCF is used in the element derivation. 

Sections 3 introduces basic terminology while Sections 4-6 describe the CCF in general 
terms. The formulation is then applied to the three-dimensional beam element in Sections 
7—11. Section 12 present numerical examples. 

3. NONLINEAR MATRIX EQUATIONS 

In this section we summarize the discrete governing equations of a geometrically nonlinear 
structure expressed in terms of a set of generalized coordinates q that for the moment are 
left unspecified. The resulting quadratic forms in q contain deformation-dependent kernel 
matrices collectively called stiffness matrices. This deformation dependency changes with the 
variational level. In the sequel we examine variational levels 0, 1 and 2, otherwise identified 
as the energy, residual- force equilibrium, and incremental levels, respectively. 

Variational Level 0: Potential Energy. The internal energy U is a nonlinear function of the 
generalized coordinates formally expressable as 

U= iq^q + qV- (1) 

The component of U that is linear in q is the prestress force vector p°. The component of U 
that is quadratic and higher in the freedoms q is assigned the kernel K. U . This is a symmetric 
matrix with dimensions of stiffness, called the energy stiffness. The total potential energy is 

J = U- q T p, (2) 

where p is the vector of applied generalized forces conjugate to q. Throughout the present 
work the applied loads p are assumed to be conservative and deformation independent. 

Variational Level 1: Residual Force Equilibrium. The first variation 8U = f T Sq of the strain 
energy defines the internal force vector f = dU/d q. Under certain conditions studied later 
this vector may be expressed as 
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This relation defines the secant stiffness matrix K r , which (if it exists) is generally unsym- 
metric. The force residual is the difference between internal and external forces: 

r=f-p. (4) 

Setting r to zero gives the discrete equilibrium equations. 

Variational Level 2: Incremental Equilibrium. The force equilibrium equation r = 0 is 
nonlinear in q. This equation is usually treated numerically by continuation procedures that 
search for solutions in the neighborhood of a previously computed equilibrium point. To 
implement this technique the residual is expressed as a function of q and a continuation 
parameter A that parametrizes the applied forces: 

r(q,A) = f(q)-p(A) = 0. (5) 

Differentiating with respect to A yields the first-order rate equations 

r' = Kq' - p' = 0, (6) 

where primes denote differentiation with respect to A. Multiplying by d\ and converting the 
d' s to A’s gives the popular incremental form K Aq = Ap. 

The tangent stiffness is fundamental in incremental-iterative solution methods and stabil- 
ity analysis, whereas the secant stiffness (by itself or embedded in the internal force vector f) 
is important in pseudo-force methods. The energy stiffness enjoys limited application per se 
but has theoretical importance as source for the other two. In the sequel we use the notation 
K level to collectively designate these three matrices. 

4. CORE PHASE OF CCF 

The core phase of the CCF establishes nonlinear response equations at the particle level, 
using the displacement gradients as degrees of freedom. The resulting equations depend 
on the mathematical model under consideration — bar, beam, plate, shell, 3D continuum, 
etc. — insofar as the form of the internal energy density, but are otherwise independent of 
finite element discretization decisions. 

Under the effect of conservative loads the structure displaces from a reference configuration 
Co, with particle coordinates JV, , to a variable current configuration C, with corresponding 
particle coordinates X { . The particle displacements u; = x, - — AT,- are collected in u. Let the 
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state of strain at C be characterize by n„ strain components e, collected in array e, and let 
the n, conjugate stresses be Si, collected in array s. 

We assume that strains stay small so that the structure remains linearly elastic. Using 
the summation convention the elastic stress-strain relations may be written 

Si = s'- + Eije-i, or s = s° + Ee, (7) 

where s® are stresses in Co (stresses that remain if — 0, also called prestresses) and Ejj = Eji 
are elastic moduli arranged as a symmetric matrix E in the usual manner. The strain energy 
density may be written 

U = e,s° + ^eiEijej = e T s° -j- ie r Ee. (8) 

The total strain energy U is obtained by integrating (8) over the structure volume: U = 
fv 0 H the integration taking place — as can be expected in a TL description — over the 
reference configuration geometry. 

Introduce now the n g displacement gradients g mn = du m /dX n . These are alternatively 
identified as <7, (t = 1,2, ...n ff ) so they can be conveniently arranged in a one-dimensional 
array g. Following Rajasekaran and Murray 2 and Felippa 4 , assume that the strains e* are 
linked to the displacement gradients through matrix relations of the form 

e « = h fg + |g T H,g, i = l,2,...n, (9) 

where h, and H, axe arrays of dimension n g x 1 and n g x n g , respectively, with H, symmetric. 
If the Green-Lagrange (GL) strain measure is chosen, all H,- are independent of g, a restriction 
enforced throughout this work. 

Denote by S u , S r and S the energy, secant and principal tangent core stiffness matrices, 
respectively, where the qualifier “principal” is explained in Sections 5-6. Symbols $ and \E , ° 
denote the core counterpart of the force vectors f and p°. With this notation the first and 
second variations of the strain energy density can be expressed as 

SU = Sg T (S u g + 9°) + ig T SS u g = Sg T (S r g + *°) = Sg T &, (10) 

S 2 U — (g T S r Sg + Sg T SS r g + (S 2 g) T & = Ag r S Sg + (S' 2 g) T &. 

These variational equations implicitly determine S r , $ and S from S u and S^ 0 . 


( 11 ) 
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5. SPECTRAL FORMS 

General expressions for S u , S r , S and $ are given in Felippa and Crivelli 1 where it is noted 
that both S u and S r may have arbitrary coefficients. More compact expressions called spectral 
forms (because of their formal similarity with the spectral decomposition of a matrix as the 
sum of rank-one matrices) can be presented at the cost of some generality. Introduce the 
vectors 

Ci = ^ + iH.g, b< = ^ = h, + H,g. (12) 

Then e,- = cfg, s.-e,- = Eijg T cf c ; g, and the spectral forms, with the summation convention 
implied throughout, are 


s" = Eij dCj + s“H„ 

(13) 

S r = £ 0 c i c, + l(4+s i )H i , 

(14) 

$ = Sjb„ 

(15) 

d<& 

~q— = Eijbibj + SjHj = Sm + Scrp. 

(16) 


The decomposition (16) of the principal tangent core stiffness S should be noted: 

1. The material stiffness matrix S m = (dsi/dg)bi = EijbibJ , which depends on the 
constitutive coefficients and displacement gradients. 

2. The principal geometric stiffness matrix Sop = Sj(3 bj/dg) = s,H,, which depends on 
the current PK2 stresses. 

The second variation (11) of the internal energy density has S = Sm + S gp as kernel 
of the quadratic form in 8 g. The core internal force $ also appears in the inner product 
(6 2 g) 3>. This second term may either survive or drop out depending on the relation of g 
with the target physical or generalized coordinates chosen in the CCF transformation phase, 
as discussed in the following section. 

If the term drops out, S becomes the tangent core stiffness and the qualifier “principal” 
becomes superfluous. If the term survives, it contributes to what we call the complementary 
geometric stiffness. 

6. TRANSFORMATION PHASE OF CCF 

The core stiffness matrices and force vectors given in (13)— (16) pertain to each material 
particle of the structure. They can be used to construct physical stiffness matrices and 
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force vectors through a process involving core-to-physical transformation and integration. 
This transformation phase of the CCF converts the core equations to discrete equations in 
terms of physical or generalized degrees of freedom collected in vector q, which defines the 
mechanical behavior of a reference volume Vo (which may specialize to a line or area domain). 

The transformation process yields stiffness matrices K^*'**, internal force vectors f and 
prestress forces p°. The effect of these transformations depend on the nature of the relations 
between the source core variables g, and the target variables q, which for the moment are 
assumed to be independent. Three possibilities exist. 

Case I. Linear relation between g and q. We have g = Tq, or in indicial form y, = Tijqj, 
where T is a transformation matrix independent of q. Invariance of U and its two variations 
taken over the reference volume Vo yields 

K Uvel = f T T S Uve! T dV Q , f= / T T *dV 0 , p°= f T T *° dV 0 . (17) 

Jv 0 J Vi, JVq 

The stiffness transformations at all levels are congruential, which gives the standard CCF 
its name. This is the case for continuum and structural finite elements that possess only 
translational degrees of freedom. Two examples of such elements are given in the survey by 
Felippa and Crivelli 1 . 

Case II. Nonlinear algebraic relation between g and q. We have g = g(q) or in index notation, 
9i = 9i(9j)- Differentiating with respect to the g, variables yields 

8gi = ^ Sqj = TijSqj, or <5g = T<5q, 

L Q 0 ( 18 ) 

8 2 9 i = Q q ^Q qk % S 9k + 0 — F ijk Sqj 8q k , or S 2 g = (F £q) 8q, 

where (F£q) is the matrix Fiji fe 8q k — Fij k Sqj; F being a cubic array. The second term in the 
expansion of S 2 qi vanishes because the qi are assumed to be independent variables. Invariance 
of 8 2 U yields the tangent stiffness transformation 

K = J y {t T (S m + S gp )T + Q} dV 0 = K M + K GP + K GC = K M + K G , (19) 

where the entries of Q are Qj k = Q k j = Fij k $i; note that Q is symmetric because F,j k = Fi k j. 
Integration of Q over Vo yields the complementary portion K GG of the geometric stiffness 
K g . The force vectors f and p° are given by the last two expressions in (17) with T defined 
in (18). 
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What happens to ht u and K r ? They can be obtained, somewhat artificially, by defining 
g = Gq where Gij = gi/qj with some care given to the limits qj —* 0. Then K u = 
f v<) G r S U G dV 0 and K r = f Vo T T S r G dVo- Because in general T ^ G one cannot expect 
symmetry in the secant stiffness. 

This case occurs in finite elements with degrees of freedom that are fixed-axis rotations, 
because plane rotations are integrable. Examples are provided by two-dimensional beams, 
and plane stress elements with drilling freedoms with only in-plane motions considered. 

Case III. Differential relation between 8 g and £q. In this final case only the variations of g 
and q can be connected: 


% = Tij 8qj, or 8g = T<Sq, 

6 2 g, = (qj 6qt = F, lt 6 qj Sq t or S 2 g = (F6q) «q 


( 20 ) 


Transformation equations (19) still applies for K and (17) for f and p°. But no integral 
g = g(q) 85 in Case II exists. Consequently K. u and K r , which need a “secant” relation 
g = Gq, cannot be constructed. Furthermore Q is not necessarily symmetric; a condition 
for that being Fijk = Fikj or dTij/dqu = dTik/dqj. Case III occurs when three-dimensional 
finite rotations appear as degrees of freedom, as in the present development. 

Up to this point the q have been assumed to be independent variables. But for complicated 
elements, such as the present one, the CCF transformations are more conveniently applied 
in stages because of the reasons noted in the Introduction. The target variables in one stage 
become the source variables for the next one. 

What happens if the q are intermediate variables in a transformation chain? If q are linear 
in the final independent degrees of freedom v, all previous formulas hold because Case I applies 
to the remaining transformations, which are strictly congruential. But if the q are nonlinear 
in v or only a nonintegrable differential relation exists, term {dgi/dqj) 6 2 qj = Tij S 2 qj in the 
second of (18) survives. The net effect is that the complementary geometric stiffness acquires 
a higher order component, implicitly defined as the kernel of 


X 0 


QiTij 6 2 qj dV 0 , 


( 21 ) 


This term cannot be resolved (meaning the explicit extraction of its stiffness kernel) until 
the transformation chain reaches downstream variables that either are the final degrees of 
freedom (and thus independent), or depend linearly on such. It is difficult to state detailed 
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rules that encompass all possible situations. Instead the treatment of the 3D beam element 
transformations in Sections 9—11 illustrates the basic techniques for “carrying forward” terms 
such as (21). 

7. KINEMATICS OF BEAM DEFORMATION 

The beam kinematics is based on two assumptions: the TL description, and the Timoshenko 
beam hypothesis: plane sections remain plane after deformation although not necessarily 
normal to the deflected longitudinal axis. The beam is isotropically elastic with Young’s 
modulus E and shear modulus G. The reference configuration of the beam is straight and 
prismatic although not necessarily stress free. A local reference frame n, is attached to it, 
with ni directed along the longitudinal axis (the locus of cross section centroids). Axes n 2 
and n 3 are in the plane of the left-end cross section; these will be eventually aligned with the 
principal inertia axes to simplify some algebraic expressions. Along these axes we attach the 
coordinate system (Yi, X?, X 3 }. This description is schematically shown in Figure 1. 

We further define a set of moving frames, denoted by {a!,a 2 ,a 3 }, parametrized by the 
longitudinal coordinate X\, Initially these frames coincide with {nj, n 2 , n 3 }, and displace 
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rigidly attached to the cross-sections of the moving current configuration, as depicted in the 
same figure. 

A beam particle originally at , -X 2 , X 3 ) displaces to 


x(x) = x„(x,) + R I ’(jr 1 )c(Jr 2 ,x 3 ), 


( 22 ) 


where x<j describes the position of the centroid of the given cross-section, R is a 3-by-3 
orthogonal matrix function that orients the displaced cross section, and £ is the cross-section 
coordinate system < T = [ 0 X 2 X 3 ]. The displacement field is 


u = x- X = Uo + (R T - I)C 


(23) 


where u 0 (JYi) = Xo(Xj ) — X(Xj) is the centroidal displacement (see Figure 1). 

In what follows 3x3 skew-symmetric matrices are consistently denoted by placing a tilde 
over the symbol of their axial 3-vector symbol; for example: 


a = spin (a) = 


0 

03 

-a 2 ' 


f ai l 

~ a 3 

0 

a i 

, a = i 

02 

a 2 

-a 1 

0 


l 03 , 


> = axial (a). 


(24) 


The skew-symmetric curvature matrix k is defined by k = R(dR T /cLY 1 ), which is the rate 
of change of the orthogonal rotation matrix R with respect to the longitudinal coordinate. 
The curvature vector is k = axial (k). We shall also require later the variation of angular 
orientation 8 0, defined as the axial vector of the skew matrix R£R^: 


<50 = R 6R t = -<5R R t , <50 = axial (<50), 


(25) 


To define the finite strain measures for the beam, we shall use the displacement gradient 
matrix 



’9 11 

512 

5l3 

G = 

9 21 

9 22 

9 23 


.531 

9 32 

9 33. 


[ 8i 62 83 1 QX . ^ 


(26) 


where I is the 3-by-3 identity matrix, and g, are 3-component gradient vectors defined for 
convenience. The 9-component gradient vector is g T = (gj gj g[). Next introduce the 
three unit 3-uples hy, j = 1,2,3, such that the j th term of hy is equal to one and the other 
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two terms are equal to zero. With the help of these quantities, explicit expressions for the 
displacement gradient vectors g can be given as 


g, = ^ + R ^<=^ + R^ K , 

g 2 = (R T - I)h 2 , g 3 = (R r - Ijh,. 


( 27 ) 


The GL finite strain measure is defined as ^ (G + G^ + G^G). The only nonzero components 
of this tensor are 

*n= hf gl + igfH gl , 

712 = 2ei 2 = hJ’gj + hfg 2 + gjHg 2 , ( 28 ) 

713 = 2 e 13 = hJ’gj + hfg 3 + gfHg 3 , 

where H is here the 3 x 3 identity matrix. By appropriate matrix expansion these strains 
could be expressed in a form befitting the general expression (9) that involves the 9- vector g. 
Observe that the orthogonality of R gives 

e 22 = hjg 2 + \Z2Z2 = R22 ~ 1 + 1(^21 + (R22 ~ l ) 2 + R% 3 ) = R22 ~ 1 + |(2 - 2 i? 22 ) = 0 , 

2e23 = h 2 g 3 + hjg 2 + g]T g 3 = R32 + -^23 + -^ 21-^31 + R22R32 ~ R32 + R23R33 — R23 = 0 ) 

(29) 

and similarly e 33 = 0 . This confirms that the only nonzero strains are ( 28 ). 

The nonzero strains may be rewritten in a more physically suggestive form: 


eu = eft + e/, 




7 = 712 + 7 . 3 , ^ = R ( h *+^), 

712 = 72 + 72 = hJV + h 2C T/ «, ^ 

713 = 73+73 = h 3<l> + h%C T K - 


Here ej, e f are stretching and flexural normal strains, 72 an d 73 represent bending-induced 
shear strains, 7 2 , 7 3 are torsion-induced shear strains, <f> is the angular distortion vector, and 
K e is the effective bending curvature defined as K e = <f>K. The last term in e/ represents a 
squared-curvature contribution to flexure, which can usually be neglected. 

From ( 8 ) and the fact that e 2 2 = e 33 = e 23 = 0 , the strain energy stored in the current 
configuration is 

Uss if r f UdA ° dXu with W=i[£;e? 1 +G( 7 1 2 2 +7 1 2 3 )]-f-s; i e 1 i+s? 2 7 12 + s; 3 7 1 3. 

J Lo J Aq 

( 31 ) 
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8. CORE BEAM EQUATIONS 


To link up with the general formulation of Sections 4-5, we arrange the nonzero GL strains 
(28) as ej = cn, ti = 712 and 63 = 713, and the associated PK2 stresses as sj = su, 
S2 = S12 and 53 = 513. The elastic moduli axe E\\ = E, E22 — £33 = G, others zero. The 
constitutive equations are sn = + Ee u , si 2 = $12 + ^712, and 513 = 5° 3 +G713. Although 

the displacement gradient vector g has 9 components, the formulas that follow axe usually 
expressed in more natural form by using the 3-component subvectors g- defined in (27). 

The spectral core stiffnesses can be compactly expressed in terms of the vectors c* = 
hj + ^Hgi and bj = h,- + Hgi for i = 1, 2, 3, where no subscript is needed in H = I. Applying 
(13) we obtain for the spectral energy stiffness 


S u = 


[ESY + G(Sl + S%) 
G Sj 7 * 

G Sj 7 

G Sj 7 - 


[5?iH 

-?2« 

-S»H‘ 

G sY 

0 

+ 

5? 2 H 

0 

0 

Gsf 

0 

GSf. 



0 

0 


(32) 


where S^ 7 = CiCj\ S^ 7 = C2C 2 ', S^ 7 = c$cj, S^ 7 = c 2 Cj' and S^ 7 = c 3 c^\ At the force residual 
level we obtain for S r a form similar to (32) except that the prestresses s®j , j = 1, 2, 3 have 
to be replaced by the midpoint stresses ^(s^ + siy) in the second matrix. 

The internal force vector conjugate to 6 g is $ = S r g + $° = + 4> r , in which 



( Si 2 b 2 + 5i3b3 
= \ 3l2bi 

5i3bi 


(33) 


represent the contribution of the normal and shear stresses, respectively. 

The principal core tangent stiffness matrix S = Sm + S gp is obtained from (16). The 


material stiffness is 


Sm = 


ES 1 +G(S 2 +S 3 ) 

gs 4 

gs 5 ‘ 

GSj 1 

GS l 

0 

G S 5 r 

0 

CD 

l—* 

1 

S3 = b 3 bj\ S 4 = 

b 2 bf 

and S 


(34) 


geometric stiffness is 



s n H 

5 12 H 

si 3 H‘ 


■(*?, +Ee„)H 

(S?2+G 7 12)H 

(5?2+G 7a3 )H' 

S GP = 

5 12 H 

0 

0 

= 

K, + gtizJh 

0 

0 


,s 13 H 

0 

0 


. ( 3 i 2 + ^7i3)H 

0 

0 


( 35 ) 
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rp 

The contribution of (<5 2 g) $ to the complementary geometric stiffness depends on the 

target variables in the ensuing transformation phase described in the next three sections. 
This phase is carried out in three stages: 

1. From particle displacement gradients to generalized gradients w at each cross section. 
An integration over the cross section area is involved. 

2. From generalized gradients w to cross-section orientation coordinates z. The rotational 
parametrization is introduced at this stage. 

3. From cross-section orientation coordinates to finite-element nodal degrees of freedom v. 
An integration over the element length, as defined by the displacement interpolation, is 
involved. 

The first two transformation stages are summarized in Tables 1 and 2, which together 
also serve to define notation. 

9. TRANSFORMATION TO GENERALIZED GRADIENTS 

The first set of target variables are the generalized gradients w(AT] ) at each reference cross 
section defined by the longitudinal coordinate X\ . The components of w are indirectly given 
through their first variation: 


* w = {3xf axr ®} r - (36) 

where £0, defined in (25) measures the variation of angular orientation. Because this is not 
generally integrable for three-dimensional motions, it is not possible to express © as a unique 
function of the displacements. 

The variation of gj is 


= 


d<5u 0 


t r 


(37) 


where we used the relation 2 ® 8k — d8Q J dX\ -f- kSQ . On using the commutative law ab — b a 
and Jacobi’s identity ab = ab - ba we may rewrite (37) as 


<*gi = 


d6u 0 T ~rdS& 


dXi 


+ R C 


-I- R t k t C6Q. 


dXi 


(38) 


For the other gradient vectors we have <5g 2 = <5R T h 2 = R 7 ’^0h 2 = R T h^<50 and 8g 3 = 
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GO 
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R h 3 6®. These relations can be collected in matrix form as 


( s Si 


*g = < 


<5g 2 




I R t C T R t k t C 
0 0 R T h^ 

0 0 R r h^ 


( dS Uo 1 


wr 



dS® 

ix; 

► = 

w 2 

6® 

i > 


w 3 . 


Sw — W <$w, (39) 


where I is the 3-by-3 identity matrix and W, axe 3-by-9 matrices. The second variation of g, 
which is required for the complementary stiffness, is 


s 2 g, = R T 6Bi T ~ + R T (efiis@ + ^. T fsen 

u X i dX i 

+ 6@ t Z((0R + 6 2 R T ( T K + R T ( T S 2 K, <40) 

S 2 g 2 = 6 2 R T i 2 , 6 2 g, = 6 2 R r i 3 


At this point it is appropriate to introduce the following section resultants: 

r = A<r b + v°, 

Q = 0 iAt+Q°, 

M ff = EI s K e + Ml, 

At r — faGIpK + Al^, 

Here V, Q, and Ai r are axial forces, transverse shear forces, bending moments and 
torsional moments, respectively, at the current configuration C\ P°, Q° , and Eire 

similar quantities at the reference configuration Co; /?i and axe transverse- shear and torsion 

coefficients, respectively, that compensate for the actual shear stress distributions; and I 5 and 
Ip are the cartesian and polar inertia tensors, respectively, of the cross section. Should the 
axes X 2 and X 3 be aligned with the principal inertia axes the latter simplify to 



"0 

0 

0 ' 


I22 + 133 

0 

o' 

1 s = 

0 

h 2 

0 

, I p = 

0 

0 

0 


0 

0 

1 

<rs 

•4? 


0 

0 

0 


s b = Ee b , 

T=T2 + T3, T2 = G72II2, T 3 = C? 73 h 3 , 


-L 


CC dA, K t — 


(41) 


~T ~ ~ T 

Ip = h 2 Is h 2 + h 3 Is h 3 . 


Because the relation between g and w is of differential type the applicable transformation 
rules are those of Case III in Section 6, and no energy or secant stiffness survives. Thus only 
the internal force vector H and tangent stiffness S associated with w axe derived below. 



Core Level Section Gradients Section Orientation Physical DOF 
Particle Section Corss-Section Whole Element 


20 
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Internal Force Vector 

The generalized internal force vector is 


n = f A wT * dA = Yl J siW T bidA = n <r + n r , 


(43) 


where 71 a and 7t r are the contributions of the normal and shear stresses respectively. Detailed 
calculations 26 result in the following expressions: 


(44) 


Observe that the term VH T <f> corresponds to the internal force obtained for a Total La- 
grangian truss element. 4 In (44) and in the sequel we neglect the squared-curvature contri- 
bution to flexure, consistent with the approximation introduced in equation (30). 

For small deformations we have R « I, 0 « hj, K e as k and kM a « 0. If these 
approximations are made, 



'R T {V<j> + kM 9 y 


R t Q 


* T M„ 

i 7£ r — 

M r 




<f> Q + k T Mr 


- V h, ■ 


- Q ' 

~T 


Mr 

h, M a 

? — 

0 


~T 


h i Q 


7t„ = 


These expressions resemble the classic linearized theory equations. 26 

Tangent Stiffness 

For the tangent stiffness we have the decomposition introduced in Section 6: 


(45) 


S — Sm + Sgp + &GC- (46) 

Furthermore, since w is nonlinear in downstream variables, the complementary geometric 
stiffness splits into two components: 

$gc = Sgcw + &gcz, (47) 

where S GCw and S C Cz contains terms that depend on the first and second variations, respec- 
tively, of R and k. The notation is suggested by the fact that S G Cw can be merged into S GP 
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to yield the geometric stiffness S G w = $GP + $GCw, which is associated with the generalized 
gradients w and independent of the rotational parametrization selected in the next set of 
target variables z. On the other hand, the kernel S GGz cannot be extracted at the w level 
and must be earned forward to the z level because it is parametrization dependent. 

Each of the components in (46)-(47) may be expressed as the sum of two contributions, 
one from the normal stresses and one from the shear stresses: 

Sm = SMa + $Mt, SgP — SgPo + SgPt, $GCx = SgCxit + $GCtt, x = w,z. ( 48 ) 
Material Stiffness 

The generalized core material stiffness is given by the congruential transformation 

SM = Sa WTSmW<1A = YjJ A Ei W T b,bf W dA = Sm<t + S Mt . 

Carrying out the algebraic manipulations one obtains 26 

R T + k t Isk)R R t kI s ^ R t »cIsK, 


( 49 ) 


= E 


symm 


<f> I s<f> <j> Is^e 

Kl I S X e 


( 50 ) 



' u4R t Ij_R 0 

ar t i ± 4> 








'0 

0 

O' 

Smt = otG 

Ip 

Ipk 

, in which lx = 

0 

1 

0 


symm 

A<j> I '\i_<j) + k T Ipk 


0 

0 

1 


( 51 ) 


Observe that the contribution R? <j)(f) T R is the core material stiffness of a Total-Lagrangian 
bar element. 4 

Geometric Stiffness due to Normal Stresses 

It is convenient to work out together all geometric stiffness terms produced by the normal 
stresses, i.e. 

&G<t = &GP<t + &GCw<t + &GCZO = &Gw<r + &GCz(r- ( 52 ) 

The appropriate definitions are 

Sgp<t = / sji W^HWj dA, 

Ja 

Sgc<t= [ s n b l S 2 gdA = Sw T S G cw < TSw + ^(S 2 R,S 2 K), 

Ja 


( 53 ) 
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where T contains Sgcz as z level kernel. Carrying out the algebraic manipulations one 
obtains 26 


$Gwa — SgPo + ^GCwtr = 




VI R M 


0 


symm 


M<,4> 

<$>M.„k + k T 


(54) 


The term VI corresponds to the well known generalized core geometric stiffness of the 
Total Lagrangian bar element. 4 

The higher order term in (53) may be expressed as 

V t ,(6 2 R,8 2 K) = Ml4>S 2 K + <f> T RS 2 R T kM IT = Sz T (v^M^ + VikM^^ 8z. (55) 
Thus 

Sgczo = V(0 r M<r) + U^Alo-; (f>) (56) 

Since the next-level target variables z include the finite rotation parametrization, matrices V 
and U depend on that choice. They are the source of unsymmetries in the stiffness matrices 
when certain rotational parametrizations are adopted, such as the incremental rotation vector. 
For the rotational vector defined in (63), these matrices are symmetric. 

Geometric Stiffness due to Shear Stresses 

The contribution of the shear stresses to the geometric stiffness is 


$Gr = «5cPr + &GC wr + «5gc zr — WT + Sgc zr ■ 


(57) 


The appropriate definitions are 

Sgpt = j si 2 (W?HW 2 + W 2 HW,) + si 3 (W^HW 3 + W 3 HWi)dA 
Sgct = f (si2b 2 +si 3 b 3 )6 2 g dA = 8w t Sgcwt^ + V t (8 2 R,6 2 k). 

J A 

Carrying out manipulations one obtains the surprisingly simple form for Sa wr 

0 


$Gwt = SgPt + SgCwt = 


0 R t Q T ' 
0 0 

L symm 0 


(58) 


(59) 
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The terms due to the second variation of g become 

Tr = Q t S 2 R$ + M?6 2 k. 
The stiffness kernel carried forward to the z level is 

Sgczt = V(M r ) + U(Q;$). 


10. TRANSFORMATION TO THE ROTATIONAL VECTOR 

The next transformation stage passes from w to z, which is a 9- vector of generalized displace- 
ments, also associated with a beam section, which embodies the parametrization of the cross 
section rotation: 

/ du 0 da \ T ( d6 u 0 dSa \ T 

Z ~{dX 1 dX x a j ’ \dX 1 dX x 6a ) ’ (62) 

Here a denotes the rotational vector parametrization defined by the standard formulas 

a = axial a, R = exp(d T ), (63) 

and which may be extracted from R by 

a = log R = — axial (R r — R), r = i|| axial (R r — R)||. (64) 

Because only the variations of w are known the relation between w and z is also of 
differential type: 


I 0 


n v , x dY(a) 

6w = Z6z, or iw= 0 y O) ~3XT S igp. ► . 

0 0 Y(a) 


d6uo 

1XT 


in which 


Y(a) = 


sin a 


( ^ sin |a N 

aa T 

1 — cos |a| 

V |a| , 

) |a 2 

aj 2 


Applying the transformations (65) we find for the internal force and the material and 
principal-geometric components of the tangent stiffness matrix: 

f, = Z t (R„ + Hr), Km, = Z r (S„)Z, Kgp, = Z T (5 G „)Z. 


(67) 
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The materialization of the geometric stiffness terms Scczv and $GCzr for the rotational 
vector is more complex. We state here only the final result: for F T 



1 

O 

O 

0 

1 


1 

O 

O 

O 
1 

U(Q;*) = 

0 0 

symm U r 

, V(M r ) = 

0 V[ 
symm V 2 


(68) 


where 


U r = Cl Q r $I + c 2 (Q$ r + $ t q) + c 3 (Q T *aa T + Q T a$I + Q T $a T + a$ r Q) 
c 5 (Qoc T 4- aQ 7 ^ + at r Q (a& T + + Q r ara r <&I^ + 

6c<&ctac T + C6 Q T aa 7 $aa r 

T --w- 

V[ = c 2 M T + c^aM^ + c 5 aa T M r + c 7 {M T ot T + a T M T l) + c s a T M T aa 7 


+ 


V 2 r = -C 3 



C4 


da 

lx 7 


M T aa T - 1- 


da 


da 


M T a l + aM x r — + a 1 — -M r I + 


dX x 


dX x 


M r aa + C7 


da 

lx T 


Mi + Mr 


da 

lx 7 


+ 


dot ( 7 m rr* \ dot rn r 

~TTr~ ar ( aM T + M r ot + a M T 1) c^-r—- aa M r ota J 

dX 1 V / dX\ 


in which 


Ci = 


sin a 


a 


c 2 = 


ci -I- 3 c 3 

C4 — 9 , C5 


a" 


c 7 = 


1 + Cl 


c 8 = 


1 — cos a; 

“ ’ 

Ci + 2 c 2 
^ : 
3c 3 — 2c 2 


c 3 = 


sin a — a cos a 


C6 = 


c 9 


ar 

C 3 + 4c 5 

a 2 

C5 — 5c 8 


(69) 

(70) 


(71) 


(72) 


a* a* 

A similar approach can be taken with (56), which defines T„. The tangent stiffness matrix 
can be obtained by superposing all contributions. 


11. TRANSFORMATION TO FINITE ELEMENT FREEDOMS 

The final stage introduces a finite element representation for the degrees of freedom. The beam 
or beam assembly is divided into a set of two-node finite elements. Each of these nodes has 
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three displacement degrees of freedom and three rotational degrees of freedom corresponding 
to the three (Jfi, X2, X3) components of the rotational vector cc. Each element in turn has 
twelve freedoms which are collected in the array v T = {d n a n } T where d n collects the six 
translational freedoms while a„ collects the six rotations. The cross-section state vector z is 
approximated inside each element by 



(73) 


where N is a matrix of linear shape functions. Since Sz = D 6v the final internal force vector 
f and tangent stiffness matrix K of each element are obtained as 

f = J D r f z dL, K = J D t K zDdL. (74) 

where L is the element length. 

The choice of shape functions for the rotational vector poses some subtle questions. In 
small-deflection analysis it is common practice to select all Timoshenko beam shape func- 
tions to be linear in Xi . This choice obviously enforces nodal compatibility while preserving 
constant curvature states. But for finite deflections a linear interpolation for the rotational 
vector components cannot exactly represent a constant curvature state unless the rotations 
are about a single axis (plane rotations). The same is true if the rotation matrix R(ATi ) is in- 
terpolated linearly. On the other hand, linear interpolation of Euler parameters does preserve 
the constant curvature state. This motivated a development of an interpolation scheme that 
starts from the 4 Euler parameters e,(Xi), i = 0, 1,2,3, = 1 that orient the normal of 

a cross section at X\. These are collected in the 4- vector £ = { e 0 <?i £2 £3 } T . Given the 

eight end values e(0) and c(L) the interpolation that can copy a constant curvature vector k 
is found to be 26 

«>— ®(-5S <”> 

where £ = ^kXi, Cl = k — \/kTk. The constant curvature vector can be extracted 

from the end values through the formula 

K = M [0® - 2'o(L)l) e(0) - (e(0) - 2e„(0)l) t(0)] , (76) 

This interpolation is then transformed to the variations in terms of the rotational vector. 
Details are provided in Crivelli’s thesis 26 . 
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12. NUMERICAL EXAMPLES 
The proposed formulation is applied to a set of six problems 
Small Deflection Analysis of a Cantilever Beam 

This example evaluates the behavior of our formulation for small rotations and displacements, 
in which case it should recover the full linear response. A two-dimensional cantilever beam 
has length L = 5, cross-section area A = 4.8 x 10 -3 and moment of inertial I = 4.45 x 10 -5 , 
with elastic modulus E = 2.1 x 10 n and shear modulus G = 8.0775 x 10 10 . The beam is 
loaded on its free end with a vertical force P = 600. Two sets of results are reported. The first 
with the formulation as is, whereas in the second we introduce the residual bending flexibility 
correction of McNeal 27 . The purpose of this correction is to accelerate the convergence 
of the linear C° bending elements, so that linear interpolation functions approximate the 
behavior of Hermitian cubics. This behavior is desirable since this thin beam can be accurately 
represented by a Bernoulli beam, for which the analytical solution is indeed cubic. 

The numeric results are shown in Table 3 for meshes with increasing number of elements. 
The second and third columns correspond to the tip deflection and rotation when no correction 
is introduced while the fourth and fifth columns are the values obtained when the residual 
bending flexibility is taken into account. 


Number of 

Tip Deflection 

Tip Rotation 

Tip Deflection 

Tip Rotation 

Elements 

w.o. correction 

w.o. correction 

w. correction 

w. correction 

1 

-2.016 X 10“ 3 

-8.026 x 10~ 4 

-2.684 x 10" 3 

-8.026 x 10" 4 

2 

-2.517 x 10" 3 

-8.026 x 10“ 4 

-2.684 x 10" 3 

-8.026 x 10~ 4 

4 

-2.643 x 10" 3 

-8.026 x 10“ 4 

-2.684 x 10- 3 

-8.026 x 10“ 4 

8 

-2.74 x 10" 3 

-8.026 x 10“ 4 

-2.684 x 10“ 3 

-8.026 x 10" 4 

16 

-2.682 x 10" 3 

-8.026 x 10~ 4 

-2.684 x 10“ 3 

-8.026 x 10- 4 


Table 3: Comparison of results for the cantilever beam 
for the small displacement case. 


As expected for a Timoshenko beam element, 28 the displacements in the first case are 
too stiff when the mesh is coarse, whereas the tip rotation is accurate. When the residual 
bending flexibility correction is used, we obtain the exact solution for any mesh. 
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Large Displacement and Rotation Analysis 


The second example is one of the few nonlinear beam problem with a simple analytic solution. 
It consists of a straight cantilever beam loaded with a moment M at its free end, as shown 
in Figure 2. 



Length L = 1000 

Bending Stiffness El = 4 10 8 


Figure 2: Straight Cantilever Beam Loaded by a 
Tip Moment. Problem Definition. 


The exact solution is a circle of radius r = EI/M. The beam has length L = 1000 and 
bending stiffness El = 4 x 10 8 while the applied moment is M = 2n x 10 5 . The mesh 
consists of 10 equally-spaced elements and the full load is applied in four increments. Using 
full Newton, convergence is attained in an average number of four iterations per increment. 
The deflected shapes for half- load and full-load levels are shown in Figure 3. 

Cantilever Beam with Two Transversal Loads 

This example consists of a two-dimensional cantilever beam with two vertical loads, one 
applied at the free end and the other one close to mid-span, as shown in Figure 4. This 
problem has been considered by several authors and an analytic solution is available, see e.g. 
Ebner and Ucciferro. 29 

The cross-section area of this beam is A = 0.2, while the bending flexibility is El = 5 x 10 6 
and the shear modulus is G = 1.153846 x 10 7 . The full load is applied in ten equal increments. 
The structure has been discretized by eight equally spaced finite elements. 
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Figure 3: Straight Cantilever Beam Loaded by a 
Tip Moment. Deformed Shapes. 



Bending Stiffness El = 5 10 6 

Shear Modulus G =1.15 10 7 
Area A =02 


Figure 4: Straight Cantilever Beam Loaded by Two 
Transversal Loads. Problem Definition. 

Table 4 gives the tip displacements and rotation for different meshes and different for- 
mulations. The results obtained by the present model are compared to those obtained by 
Cardona 24 and to those provided by the analytic solution. Deformed configurations for se- 
lected load levels are shown in Figure 5. 

Large Displacement 3-D Analysis of 45-Degree Bend 
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Formulation Elements 

2 4 8 




Longitudinal 

-28.99 

-30.26 

-30.62 

■ 

Cardona 

Transversal 

-65.86 

-66.63 

-66.87 



Rotation 

-1.1 

-1.06 

-1.05 

w 


Longitudinal 

-28.99 

-30.26 

-30.62 


Present 

Transversal 

-65.86 

-66.63 

-66.87 



Rotation 

-1.1 

-1.06 

-1.05 



Longitudinal 

-30.75 

-30.75 

-30.75 

“ 

Analytic 

Transversal 

-66.96 

-66.96 

-66.96 



Rotation 

— 

— 

— 


Table 4: Comparison of results for the cantilever beam 
with two transversal loads. 




Figure 5: Straight Cantilever Beam Loaded by Two 
Transversal Loads. Deformed Shapes. 

This example has been studied by Bathe and Bolourchi. 30 It consists of a cantilever 45-degree 
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bend lying on a horizontal plane, subjected to an end-load normal to that plane, as shown in 
Figure 6. 



Figure 6: Curved Cantilever Bend Loaded by 
Tip Load. Problem Definition. 


The bend is an arc of a circle of radius r = 100 and the beam cross-section is a square 
with sides of unit length. The material has elastic modulus E = 10 7 and a zero Poisson ratio. 
The finite element mesh consist of 8 straight beam elements. The full load P = 600 is applied 
in 6 equal increments. 

Table 5 compares the results obtained by different formulations while the deflected shapes 
corresponding to selected load levels are presented in Figure 7. 


Formulation 


300 


Load 

450 


600 


Bathe et.al 
Simo et.al 
Cardona 
Present 


22.33, 58.84, 40.08 
22.50, 59.20, 39.50 
22.14, 58.64, 40.35 
22.31, 58.85, 40.08 


18.62, 53.32, 48.39 

18.38, 52.11, 48.59 
18.59, 53.34, 48.39 


15.79, 47.23,53.37 
15.90, 47.20,53.40 
15.55, 47.04,53.50 
15.75, 47.25, 53.37 


Table 5: Comparison of results for the 45-degree bend cantilever beam 
with a transversal tip load. 
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Figure 7: Curved Cantilever Bend Loaded by 
Tip Load. Deflected Shapes. 

This table shows that the result obtained by Cardona 24 are stiffer whereas those obtained 
by Simo and Vu-Quoc 23 axe more flexible than those obtained by Bathe et.al. 30 However, 
these latter results can be taken as reference since they have also been obtained using a more 
refined mesh. It can be seen that the results of the present formulation agree well with the 
reference solution. 

Williams’ Toggle Beam 

This problem consists of two thin plane beams rigidly jointed together and clamped at both 
ends. The load P is applied at the apex of the structure. The geometry and physical 
properties are given in Figure 8. Williams 31 solved this problem analytically by taking into 
account the effects of finite changes in geometry as well as the effects of the axial force in 
the flexural stiffness and the flexural modification of the axial stiffness. He also compared his 
results to experimental data. 

This problem was first solved numerically by Wood and Zienkiewicz 32 using five Total- 
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0.386 " 


Stretching Stiffness EA = 1A55 10 6 U> 
Bending Stiffness El -927 10 3 Win 2 


Figure 8: William’s Toggle. Problem Definition. 



Figure 9: William’s Toggle. Load-Deflection curves. 


Lagrangian paralinear isoparametric elements per member. In addition, these authors used a 
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displacement control strategy to traverse the limit point. This problem was also analyzed by 
Papadrakakis 33 using two Euler-Bernoulli beams, formulated under the assumption of large 
displacements but moderate rotations. Nee and Haidar 35 solved this problem by using an 
assumed stress formulation. We have solved this problem using five elements per member 
together with the residual bending flexibility concept 27 to approximate the thin beam behav- 
ior. Comparisons of the tip deflections versus load obtained by these different formulations 
are given in Figure 9. We observe a good agreement between our results and experimental 
measurements. 

Twelve Member Frame 

This problem consists of a three-dimensional frame made of twelve members, with half of 
them laid out as an hexagon, and the other half making up the diagonals of the hexagon. 
The load is applied on the central node. The geometry and physical properties of the frame 
are described in Figure 10. 



24" 


Polar Moment J x 
Second Moments Iy 
Young’s Modulus E 
Shear Modulus G 
Area A 


0.0331 in 4 
Iz = 0.0203 in 4 
439.8 ksi 
159 ksi 
0.494 in 2 



Figure 10: Twelve-Member Frame. Problem Definition. 


This frame is simply supported and the supports are allowed to move on the plane normal 
to the load. To remove the translational rigid body modes and make the problem determinate, 
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the central node is restricted against lateral displacement. The frame has been discretized 
using one finite element for each member of the base and two elements for each diagonal 
member. 



Vertical Displacement, inches 

Figure 11: Twelve- Member Frame. Load- Deflection Curves. 

The evolution of the deflection of the apex while the load is varied is given in Figure 11. 
An incremental strategy with step control and a hyperelliptic constraint have been used to 
traverse the limit point. The results obtained by the present formulation are compared to 
those obtained by Papadrakakis, 33 Meek and Tan 34 and Nee and Haidar. 35 It can be noticed 
that the present formulation displays a slightly stiffer behavior, which can be attributed to the 
presence of the shear stress. The extra stiffness should disappear with more refined meshes. 


13. CONCLUSIONS 

We have constructed and tested a three-dimensional Timoshenko beam element based on the 
Total Lagrangian description. The element has 12 degrees of freedom: 6 translations and 
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6 finite rotations measured by the rotational vector. This particular set of nodal freedoms 
reduces to the usual choice in small-deflection analysis. This uniformity of treatment allows 
implementation in standard finite element codes without the necessity of making special 
provisions to update the rotations. Furthermore, the formulation leads to a symmetric tangent 
stiffness matrix for arbitrary motions. This attribute is particularly valuable in stability 
analysis of complex structures, as it allows the use of symmetric eigensolvers near bifurcation 
points. 

The kinematics of the deformation is described using an inertial frame attached to the 
undeformed configuration of the beam. This kinematic description has a potential advantage 
in nonlinear dynamic calculations in that the mass matrix remains fixed, with the same 
sparsity as in small- deflection analysis. The use of Green- Lagrange strains, conjugate PK2 
stresses and the absence of hazardous a priori kinematic approximations (beyond those of the 
Timoshenko beam model) effectively filters out arbitrary rigid body motions, and allows the 
beam element to capture coupling effects between stretching, bending, torsion and transverse 
shears within the elastic response regime. These abilities augur well for its future use in highly 
flexible space structures, where the effect of those couplings can be extremely important in 
stability, dynamics and control. 

The discrete equations have been derived using the staged approach of the Core- 
Congruential Formulation (CCF). In the innermost level, core equations are obtained at 
the particle level. These physically transparent equations depend only on the form of the 
internal energy density. A chain of transformations ensues in which the core equations are 
referred to three sets of kinematic variables, two pertaining to cross sections and the third 
one to the finite element nodal degrees of freedom. The choice of finite rotation measure is 
introduced in the second stage. The choice of finite element interpolation and nodal free- 
doms is introduced in the last stage. This “nesting” offers obvious advantages in fostering 
programming modularity and maintaning flexibility as regards to decision changes. 

We believe that the main contributions of the present work to computational nonlinear 
mechanics are as follows. 

1. The development of a new symmetric formulation for the analysis of the geometrically non- 
linear response of three-dimensional beams that undergo arbitrary rotations. The Total 
Lagrangian description maintains a fixed reference configuration, which is advantageous 
in many classes of problems. No special treatment of the rotational degrees of freedoms 
is required thus simplifying the treatment of boundary conditions. The symmetry and 



freedom-choice attributes simplifies the element implementation into stiffness-based finite 
element codes. 

2. The CCF methodology allows a systematic staged development of the Total Lagrangian 
element equations that maintains physical visibility. The general core equations display 
microscopic behavior, and are gradually specialized to macroscopic behavior in the trans- 
formation phase. Behavioral approximations may be injected in the initial transformation 
stages, whereas computational decisions as regards rotational parametrization and element 
discretization may be deferred to the final transformation stages. 

Extensive numerical experiments have been performed to validate and test the present 
formulation and solution methods. These problems cover a wide range of structural behavior, 
from plane to three-dimensional structures, including snap-through and nonlinear bifurcation. 
The ability of the present formulation to deal with large three-dimensional rotations and 
displacements has been demonstrated. In addition, the beam model approaches the linear 
beam behavior when the displacements and rotations are small. For very thin beams it is 
well known that Timoshenko beam models may be stiffer than Euler- Bernouilli models. The 
residual bending flexibility correction may improve this behavior when the displacements Eire 
moderate. 

Although the solution method is not a focus of this paper, it is noted that good Newton- 
Raphson convergence rates have been attained for all the tested problems. This behavior 
validates the consistency of the residual force and tangent stiffness computations. 
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